Skip to content

ENH/FIX: Resolve LinAlgError during SVD #2838

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Dec 22, 2018
Merged

Conversation

feilong
Copy link
Contributor

@feilong feilong commented Dec 21, 2018

The default lapack_driver uses divide and conquer SVD (gesdd), which occasionally will cause LinAlgError even if the data are completely fine. Switching lapack_driver to gesvd when it fails can solve this issue.

Summary

Fixes # .

List of changes proposed in this PR (pull-request)

Acknowledgment

  • (Mandatory) I acknowledge that this contribution will be available under the Apache 2 license.

The default `lapack_driver` uses randomized SVD, which occasionally will cause `LinAlgError` even if the data are completely fine. Switching `lapack_driver` to `gesvd` will solve this issue.
@@ -1193,9 +1194,12 @@ def compute_noise_components(imgseries, mask_images, num_components,
try:
u, _, _ = np.linalg.svd(M, full_matrices=False)
Copy link
Member

@effigies effigies Dec 21, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nested try/except blocks get a little hard to reason about. How about a helper function that will do the fall-back:

# Use the numpy.linalg.svd protocol
def fallback_svd(a, full_matrices=True, compute_uv=True):
    try:
        return np.linalg.svd(a, full_matrices=full_matrices, compute_uv=compute_uv)
    except np.linalg.LinAlgError:
        pass

    from scipy.linalg import svd
    return svd(a, full_matrices=full_matrices, compute_uv=compute_uv, lapack_driver='gesvd')

And then we just replace the above with:

Suggested change
u, _, _ = np.linalg.svd(M, full_matrices=False)
u, _, _ = fallback_svd(M, full_matrices=False)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And to be clear, I'm suggesting we place the above fallback_svd function at the file level.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great idea! By putting it at the file level do you mean putting it at the top of the file?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. Anywhere, really, but not nested within a class or another function. The top is fine.

@codecov-io
Copy link

codecov-io commented Dec 21, 2018

Codecov Report

Merging #2838 into master will decrease coverage by 0.02%.
The diff coverage is 44.44%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master    #2838      +/-   ##
==========================================
- Coverage   67.47%   67.45%   -0.03%     
==========================================
  Files         341      341              
  Lines       43355    43362       +7     
  Branches     5379     5379              
==========================================
- Hits        29254    29248       -6     
- Misses      13405    13409       +4     
- Partials      696      705       +9
Flag Coverage Δ
#smoketests 50.53% <11.11%> (-0.02%) ⬇️
#unittests 64.85% <44.44%> (-0.04%) ⬇️
Impacted Files Coverage Δ
nipype/algorithms/confounds.py 66.28% <44.44%> (-0.32%) ⬇️
nipype/pipeline/plugins/legacymultiproc.py 61.5% <0%> (-4.5%) ⬇️
nipype/interfaces/dynamic_slicer.py 17.47% <0%> (ø) ⬆️
nipype/interfaces/nipy/preprocess.py 45.79% <0%> (ø) ⬆️
nipype/interfaces/io.py 53.99% <0%> (ø) ⬆️
nipype/interfaces/freesurfer/preprocess.py 66.11% <0%> (ø) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b582d53...053fd61. Read the comment docs.

@chrisgorgo
Copy link
Member

Great contribution @feilong! Thanks for getting to the bottom of this!

@@ -14,6 +14,7 @@
import nibabel as nb
import numpy as np
from numpy.polynomial import Legendre
from scipy import linalg
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unneeded import.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for pointing that out!

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM! Thanks a lot for figuring this one out. Hopefully this will take care of most of the cases we see in fMRIPrep.

@feilong
Copy link
Contributor Author

feilong commented Dec 21, 2018

Great! Thanks for making these great softwares. Looking forward to more commits!

@effigies effigies added this to the 1.1.8 milestone Dec 22, 2018
@effigies effigies merged commit dcc013b into nipy:master Dec 22, 2018
yarikoptic added a commit to yarikoptic/nipype that referenced this pull request Feb 4, 2019
1.1.8 (January 28, 2019)

  * FIX: ANTS LaplacianThickness cmdline opts fixed up (nipy#2846)
  * FIX: Resolve LinAlgError during SVD (nipy#2838)
  * ENH: Add interfaces wrapping DIPY worflows (nipy#2830)
  * ENH: Update BIDSDataGrabber for pybids 0.7 (nipy#2737)
  * ENH: Add FSL `eddy_quad` interface (nipy#2825)
  * ENH: Support tckgen -select in MRtrix3 v3+ (nipy#2823)
  * ENH: Support for BIDS event files (nipy#2845)
  * ENH: CompositeTransformUtil, new ANTs interface (nipy#2785)
  * RF: Move pytest and pytest-xdist from general requirement into tests_required (nipy#2850)
  * DOC: Add S3DataGrabber example (nipy#2849)
  * DOC: Skip conftest module in API generation (nipy#2852)
  * DOC: Hyperlink DOIs to preferred resolver (nipy#2833)
  * MAINT: Install numpy!=1.16.0 from conda in Docker (nipy#2862)
  * MAINT: Drop pytest-xdist requirement, minimum pytest version  (nipy#2856)
  * MAINT: Disable numpy 1.16.0 for Py2.7 (nipy#2855)

* tag '1.1.8': (79 commits)
  MNT: Add @feilong to .zenodo, update ordering
  MNT: Update .mailmap
  MNT: Update .zenodo ordering
  Accept invitation as Zenodo release co-author (see nipy#2864)
  MAINT: Update .mailmap
  BF: allowing bids_event_file as alternate input
  MNT: Update .zenodo ordering
  MNT: Version 1.1.8
  DOC: 1.1.8 changelog
  Update nipype/interfaces/dipy/tracks.py
  Update nipype/interfaces/dipy/reconstruction.py
  MNT: Install numpy!=1.16.0 from conda in Docker
  Add FSL auto test
  remake specs
  Update nipype/interfaces/io.py
  Remove return type named tuple
  Update nipype/info.py
  STY: Whitespace, line length
  Remove out_ prefix from EddyQuad outputs
  Apply minor edits from code review
  ...
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants